ELISA
2022-08-15
Load in libraries
library(tidyverse)
library(readxl)
library(mgcv)
library(emmeans)
library(ggeffects)
library(lme4)
library(formattable)
library(ggpubr)
library(pzfx)
library(DT)
library(cowplot)
library(gratia)
library(cenGAM)
source("scripts/utils/utils.R")
modernacolors()## <environment: R_GlobalEnv>Load data
# load data
d <- read_excel("processed_data/ELISA.xlsx")Observed data visualization
Supplementary Figure 3. Impact of mRNA-1273 Prime-boost Interval on S2-P-specific Serum Binding IgG Antibody Titers Through 24 Weeks Post-boost
#-------------------------------------------------------------------------------
# Supplementary Figure 1
#-------------------------------------------------------------------------------
# color panel
panel0 <-
c(
modernadarkblue2,
modernaorange,
modernapurple,
modernagreen,
modernablue,
modernared,
modernadarkgray
)
dosepanel <- c("#cb6939", "#8264cb")
LLOQ <- 25 # lower limit of quantificatin
# trend over time
raw1 <-
d %>%
filter(interval != "PBS") %>%
mutate(interval = factor(interval)) %>%
mutate(
.,
interval = forcats::fct_recode(
interval,
"8-week interval" = "8",
"6-week interval" = "6",
"4-week interval" = "4",
"3-week interval" = "3",
"2-week interval" = "2",
"1-week interval" = "1",
"Prime only" = "0"
)
) %>%
mutate(., interval = factor(.$interval, levels = levels(.$interval)[8:1])) %>%
mutate(
day = forcats::fct_recode(
factor(day),
"Prime only" = "56",
"1 wk post-boost" = "64",
"2 wk post-boost" = "73",
"4 wk post-boost" = "87",
"8 wk post-boost" = "115",
"12 wk post-boost" = "143",
"16 wk post-boost" = "171",
"20 wk post-boost" = "199",
"24 wk post-boost" = "226"
)
) %>%
mutate(dose = forcats::fct_recode(
factor(dose),
"mRNA-1273 1 µg" = "1",
"mRNA-1273 10 µg" = "10"
)) %>%
ggplot(aes(
x = factor(day),
y = value,
color = factor(interval)
)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(position = position_jitterdodge(jitter.width = .2)) +
geom_hline(aes(yintercept = log10(LLOQ), lty = "LLOQ")) +
scale_linetype_manual(values = 2) +
facet_wrap( ~ dose, nrow = 3) +
theme_bw() + theme(panel.grid.major = element_line(colour = "grey100")) +
labs(
x = " ",
y = expression(log[10] ~ S2P ~ IgG ~ titer),
color = "",
linetype = ""
) +
scale_color_manual(
values = c(
modernared,
modernablue,
modernagreen,
modernapurple,
modernaorange,
modernadarkblue2,
modernadarkgray
)
) +
scale_y_continuous(breaks = 1:7) +
theme(
axis.text = element_text(size = 12),
axis.title = element_text(size = 16),
strip.text = element_text(size = 14),
legend.text = element_text(size = 14),
legend.title = element_text(size = 16)
)
# animal-level line plot
raw2 <-
d %>%
filter(interval != "PBS") %>%
mutate(
interval = factor(interval),
animal_id = paste("animal", animal_id),
animal_id = factor(
animal_id,
levels = c(
"animal 1",
"animal 2",
"animal 3",
"animal 4",
"animal 5",
"animal 6",
"animal 7",
"animal 8",
"animal 9",
"animal 10"
)
)
) %>%
mutate(
.,
interval = forcats::fct_recode(
interval,
"8-week interval" = "8",
"6-week interval" = "6",
"4-week interval" = "4",
"3-week interval" = "3",
"2-week interval" = "2",
"1-week interval" = "1",
"Prime only" = "0"
)
) %>%
mutate(., interval = factor(.$interval, levels = levels(.$interval)[8:1])) %>%
mutate(dose = forcats::fct_recode(
factor(dose),
"mRNA-1273 1 µg" = "1",
"mRNA-1273 10 µg" = "10"
)) %>%
ggplot(aes(
x = day,
y = value,
color = factor(animal_id)
)) +
geom_line() +
geom_jitter(position = position_jitterdodge(jitter.width = .2)) +
geom_hline(aes(yintercept = log10(LLOQ), lty = "LLOQ")) +
scale_linetype_manual(values = 2) +
facet_grid(dose ~ interval) +
theme_bw() + theme(panel.grid.major = element_line(colour = "grey100")) +
labs(
x = " ",
y = expression(log[10] ~ S2P ~ IgG ~ titer),
color = "",
linetype = ""
) +
scale_color_manual(
values = c(
"#c25dba",
"#72b24a",
"#7561cf",
"#ce9b44",
"#807dc3",
"#7a7d37",
"#4faad9",
"#c95c3f",
"#4bac88",
"#c8577b"
)
) +
scale_y_continuous(breaks = 1:7) +
scale_x_continuous(
breaks = c(56 , 87, 115, 143, 171, 199, 227),
labels = c(
"Prime only",
"4 wk post-boost",
"8 wk post-boost",
"12 wk post-boost",
"16 wk post-boost",
"20 wk post-boost",
"24 wk post-boost"
)
) +
theme(
axis.text = element_text(size = 12),
axis.title = element_text(size = 16),
axis.text.x = element_text(
angle = 45,
hjust = 1,
vjust = 1
),
strip.text = element_text(size = 14),
legend.text = element_text(size = 14),
legend.title = element_text(size = 16)
)
cowplot::plot_grid(
plotlist = list(raw1, raw2),
labels = c("A", "B"),
nrow = 2,
ncol = 1
)
Model fitting
Model selection
We use a generalized additive model (GAM) with a 9-dimensional thin plate spline basis in day. We try 4 different GAMs with varying complexity.
Model A: a GAM censored regression model that accounts for censored values in prime-only group, dosing interval group-specific smooth nonlinear trend in days, a dose-specific smooth nonlinear trend in days, animal-specific random effect, and animal_specific effect of day
Model B: a GAM regression model that exclude the prime-only group, with dosing interval group-specific smooth nonlinear trend in days, a dose-specific smooth nonlinear trend in days, animal-specific random effect, and animal_specific effect of day
Model C: a GAM regression model that exclude the prime-only group, with animal-specific random effect, and animal_specific effect of day
Model D: a GAM regression model that exclude the prime-only group, with a single day-specific smooth, dosing interval specific random effect, and dosing interval specific effect of day
Model B is the best model among the four choices.
- Model A does not fit as good as Model B which excluded the prime-only data. Note prime-only data has a different trend as other dosing interval groups with some censored observations.
- Model C that is without smoothing splines does not fit as good as Model B that has smoothing splines.
- Model B’s residual diagnostics look slightly better than Model D.
Model D has smaller
BICbut largerAICthan Model B. Based onAICand residual diagnostics, Model B is better than Model D.
#-------------------------------------------------------------------------------
# Fit the model
#-------------------------------------------------------------------------------
dtemp <- filter(d, interval != "PBS") %>%
mutate(
dose = factor(dose),
day = as.numeric(day),
animal_id = factor(animal_id),
interval = factor(interval)
) %>%
mutate(
animal_unique_id = paste0(interval, "_", animal_id),
animal_unique_id = factor(animal_unique_id)
)
# fit a gam censored model that includes prime-only group with dosing interval group-specific smooth nonlinear trend in days, a dose-specific smooth nonlinear trend in days, animal-specific random effect, and animal_specific effect of day
gam_fit_A <-
mgcv::gam(
value ~ interval + dose + interval:day +
interval:dose + dose:day + dose:day:interval +
s(day, by = dose, bs="tp", k = 9) + s(day, by = interval,bs="tp", k = 9) +
s(animal_unique_id, bs = "re") + s(animal_unique_id, day, bs = "re"),
data = dtemp,
family = tobit1(left.threshold = 1.097),
method = 'REML'
)
# fit a gam model the same as previous model but without prime-only group (interval = "0")
d_mod <- dtemp %>%
filter(interval != "0")%>%
droplevels()
gam_fit_B<-
mgcv::gam(
value ~ interval + dose + interval:day +
interval:dose + dose:day + dose:day:interval +
s(day, by = dose, bs="tp", k = 9) + s(day, by = interval, bs="tp", k = 9) +
s(animal_unique_id, bs = "re") + s(animal_unique_id, day, bs = "re"),
data = d_mod,
method = 'REML'
)
# fit a gam model without prime-only group (interval = "0") with animal-specific random effect, and animal_specific effect of day
gam_fit_C <-
mgcv::gam(
value ~ interval + dose + day + interval:day +
interval:dose + dose:day + dose:day:interval +
s(animal_unique_id, bs = "re") + s(animal_unique_id, day, bs = "re"),
data = d_mod,
method = 'REML'
)
# fit a gam model without prime-only group (interval = "0") with a single day-specific smooth, dosing interval specific random effect, and dosing interval specific effect of day
gam_fit_D <-
mgcv::gam(
value ~ interval + dose + s(day, k = 9) + interval:day +
interval:dose + dose:day + dose:day:interval +
s(interval, bs = "re") + s(interval, day, bs = "re") ,
data = d_mod,
method = 'REML'
)GAMs
Estimated smooths visualization
plot_grid(gratia::draw(gam_fit_A), labels = c("A"), vjust = 1) 
plot_grid(gratia::draw(gam_fit_B), labels = c("B"), vjust = 1) 
plot_grid(gratia::draw(gam_fit_C), labels = c("C"), vjust = 1) 
plot_grid(gratia::draw(gam_fit_D), labels = c("D"), vjust = 1) 
Residual diagnostics
# Residual diagnostics plots
list(gam_fit_A, gam_fit_B, gam_fit_C, gam_fit_D) %>%
map(~assumption_check(.x, .x$model$value))%>%
cowplot::plot_grid(plotlist = ., nrow = 2, labels = c("A", "B", "C", "D"))
AIC and BIC comparison
AIC(gam_fit_A, gam_fit_B, gam_fit_C, gam_fit_D)%>%
rownames_to_column("Model") %>%
DT::datatable()%>%
formatSignif(columns = c( 'df', 'AIC'),
digits = 6)BIC(gam_fit_A, gam_fit_B, gam_fit_C, gam_fit_D)%>%
rownames_to_column("Model") %>%
DT::datatable() %>%
formatSignif(columns = c( 'df', 'BIC'),
digits = 6)Selected model residual diagnostics
# Residual diagnostics figure
gam_fit <- gam_fit_B
p <- assumption_check(gam_fit, gam_fit$model$value)
title <-
ggdraw() + draw_label(
"ELISA residual diagnostics when model day as a continuous variable",
x = 0.03,
hjust = 0
)
plot_grid(
title,
p,
ncol = 1,
rel_heights = c(0.1, 1),
align = "v"
)
Residual diagnostics
Figure 2. S2-P-specific Serum Binding IgG Antibody Titers
#-------------------------------------------------------------------------------
# Figure 2. Spike-specific Serum Binding IgG Antibody Titers
#-------------------------------------------------------------------------------
plot1 <- ggemmeans(gam_fit, terms = c("day", "interval", "dose"))
F2A <-
plot(plot1) +
scale_color_manual(values = panel0,
labels = c("1-week interval", "2-week interval",
"3-week interval", "4-week interval",
"6-week interval", "8-week interval")) +
scale_fill_manual(values = panel0) +
scale_x_continuous(
breaks = c(56 , 64, 73, 87, 115, 143, 171, 199, 226),
labels = c(
"Prime only",
"1 wk post-boost",
"2 wk post-boost",
"4 wk post-boost",
"8 wk post-boost",
"12 wk post-boost",
"16 wk post-boost",
"20 wk post-boost",
"24 wk post-boost"
)
) +
labs(
title = "", color = "", x = "",
y = expression(log[10] ~ S2P ~ IgG ~ titer)
) +
theme(panel.grid.major = element_line(colour = "grey100")) +
theme(
axis.text.x = element_text(
angle = 45,
hjust = 1,
vjust = 1
),
legend.position = "bottom"
) +
facet_grid(~ facet, labeller = labeller(.cols = ~ paste("mRNA-1273", .x, "µg")))# plot effect of dose treating dose as numerical variable
plot2 <-
ggemmeans(
gam_fit,
terms = c("day", "dose", "interval"),
condition = c("day", "interval")
)
F2B <-
plot(plot2) +
scale_color_manual(values = c(modernared, modernablue),
labels = c("mRNA-1273 1 µg", "mRNA-1273 2 µg")) +
scale_x_continuous(
breaks = c(56 , 64, 73, 87, 115, 143, 171, 199, 226),
labels = c(
"Prime only",
"1 wk post-boost",
"2 wk post-boost",
"4 wk post-boost",
"8 wk post-boost",
"12 wk post-boost",
"16 wk post-boost",
"20 wk post-boost",
"24 wk post-boost"
)
) +
labs(
x = "",
title = "", color = "",
y = expression(log[10] ~ S2P ~ IgG ~ titer)
) +
theme(panel.grid.major = element_line(colour = "grey100")) +
theme(axis.text.x = element_text(
angle = 45,
hjust = 1,
vjust = 1
),
legend.position = "bottom") +
facet_wrap(~ facet, labeller = labeller(.cols = ~paste(.x, "week interval")))
Figure2 <-
cowplot::plot_grid(
plotlist = list(F2A, F2B),
labels = c("A", "B"),
nrow = 2,
ncol = 1
)
Figure2
Residual diagnostics when model day as a factor
# fit the model for mean comparison
dtemp <- filter(d, interval != "PBS" & interval != "0") %>%
mutate(
dose = as.factor(dose),
animal_id = factor(animal_id),
interval = factor(interval),
day = factor(day)
) %>%
mutate(animal_id = paste0(interval, "_", animal_id)) %>%
mutate(animal_id = paste0(animal_id , "_", dose))
lm_fit <-
lmer(value ~ interval * dose * day + (1 |
animal_id) , data = dtemp)
p2 <- assumption_check(lm_fit, lm_fit@frame$value)
# Residual diagnostics figure
title <-
ggdraw() + draw_label(
"ELISA residual diagnostics when model day as a discrete variable",
x = 0.03,
hjust = 0
)
plot_grid(
title,
p2,
ncol = 1,
rel_heights = c(0.1, 1),
align = "v"
)
Supplementary Figure 1. Predicted Fold-change of Spike-specific Serum Binding IgG Antibody Titers
Each dot corresponds to the estimated fold change for the contrasted prime-boost dosing interval. Error bars are multivariate t adjusted 95% CIs on the estimated contrast.
# The results indicate that the peak immune response/ titer levels are achieved at 2 weeks post prime
grid <- ref_grid(lm_fit, cov.keep = c("day", "dose", "interval"))
emm_fit_day <- emmeans(grid, ~ day | interval + dose)
#-------------------------------------------------------------------------------
# Supplementary Figure 2. Predicted Fold-change of Spike-specific Serum Binding IgG Antibody Titers.
#-------------------------------------------------------------------------------
set.seed(101)
summary(
contrast(
emm_fit_day,
list(
"2 weeks post-boost \n vs Pre-boost" = c(-1, 0, 1, 0, 0, 0, 0, 0, 0),
"24 weeks post-boost \n vs Pre-boost" = c(-1, 0, 0, 0, 0, 0, 0, 0, 1)
)
),
infer = TRUE,
adjust = "mvt",
level = .95
) %>%
as_tibble() %>%
mutate(
.,
interval = forcats::fct_recode(
interval,
"8-week interval" = "8",
"6-week interval" = "6",
"4-week interval" = "4",
"3-week interval" = "3",
"2-week interval" = "2",
"1-week interval" = "1"
)
) %>%
mutate(., interval = factor(.$interval, levels = levels(.$interval)[8:1])) %>%
mutate(dose = forcats::fct_recode(
factor(dose),
"mRNA-1273 1 µg" = "1",
"mRNA-1273 10 µg" = "10"
)) %>%
ggplot(aes(x = contrast , y = estimate, colour = interval)) +
geom_point(position = position_dodge(width = 0.6), size = 2) +
geom_errorbar(
aes(ymin = lower.CL, ymax = upper.CL),
width = .1,
position = position_dodge(width = 0.6),
size = 1
) +
scale_color_manual(
values = c(
modernared,
modernablue,
modernagreen,
modernapurple,
modernaorange,
modernadarkblue2,
modernadarkgray
)
) +
theme_bw() + theme(panel.grid.major = element_line(colour = "grey100")) +
scale_y_continuous(breaks = log10(c(0.5, 1, 3, 10, 30, 100)), labels = c(0.5, 1, 3, 10, 30, 100)) +
geom_hline(yintercept = 0,
color = "darkcyan",
lty = "dashed") +
labs(x = "",
y = expression(Estimated ~ fold ~ change ~ of ~ S2P ~ IgG ~ titer)) +
facet_grid(~ dose) +
theme(
axis.text = element_text(size = 12),
axis.title = element_text(size = 18),
strip.text = element_text(size = 18),
legend.text = element_text(size = 16),
legend.title = element_text(size = 18),
plot.title = element_text(size = 20),
legend.position = "bottom"
)
Supplementary Table 1. Statistical comparisons of S2-P IgG antibody titer in mice between different mRNA-1273 dosing intervals.
#-------------------------------------------------------------------------------
# Supplementary Table 1
#-------------------------------------------------------------------------------
grid <- ref_grid(gam_fit, cov.keep = c("day", "dose", "interval"))
em_fit <- emmeans::emmeans(grid, ~ interval | day + dose)
set.seed(100)
sum_fit <-
summary(
pairs(em_fit , reverse = TRUE),
infer = TRUE,
adjust = "mvt",
level = .95
)
# generate tables
elisa_p_table <-
sum_fit %>%
as_tibble() %>%
mutate(day = forcats::fct_recode(
factor(day),
"0 (pre-dose 2)" = "56",
"7 (1)" = "64",
"16 (2)" = "73",
"30 (4)" = "87",
"58 (8)" = "115",
"86 (12)" = "143",
"114 (16)" = "171",
"142 (20)" = "199",
"169 (24)" = "226"
))%>%
dplyr::select(contrast, estimate, day, dose, p.value) %>%
mutate(estimate = 10 ^ estimate) %>%
mutate(contrast = str_remove_all(contrast, "interval"),
contrast = str_replace(contrast, " - ", " / ")) %>%
rename(
"Pairwise comparison" = contrast,
"Fold change" = "estimate",
"Adjusted P-value" = "p.value"
) %>%
mutate(`Adjusted P-value` = round(`Adjusted P-value`, digits = 6),
`Fold change`= round(`Fold change`, digits = 6)) %>%
filter(`Adjusted P-value` < 0.05)
elisa_p_table %>%
DT::datatable(caption = "Statistical comparisons of S2-P IgG antibody titer in mice between different mRNA-1273 dosing intervals")Supplementary Figure 4. Statistical comparisons of S2-P IgG Antibody Titers between different mRNA-1273 dosing intervals(Supplementary Table 1 visualization)
Each dot corresponds to the estimated fold change for the contrasted prime-boost dosing interval. Error bars are multivariate t adjusted 95% CIs on the estimated contrast.
sum_fit %>%
as_tibble() %>%
mutate(contrast = str_remove_all(contrast , "interval")) %>%
mutate(contrast = str_replace(contrast, " - ", " weeks / ")) %>%
mutate(contrast = ifelse(
str_detect(contrast, "1"),
paste(contrast, "week"),
paste(contrast, "weeks")
)) %>%
mutate(day = as.factor(day)) %>%
mutate(
day = forcats::fct_recode(
factor(day),
"Prime only" = "56",
"1 wk post-boost" = "64",
"2 wk post-boost" = "73",
"4 wk post-boost" = "87",
"8 wk post-boost" = "115",
"12 wk post-boost" = "143",
"16 wk post-boost" = "171",
"20 wk post-boost" = "199",
"24 wk post-boost" = "226"
)
) %>%
rename(
"Pairwise comparison" = contrast,
"Fold change" = "estimate",
"Adjusted P-value" = "p.value",
"Day" = "day"
) %>%
mutate(dose = forcats::fct_recode(
factor(dose),
"mRNA-1273 1 µg" = "1",
"mRNA-1273 10 µg" = "10"
)) %>%
ggplot(aes(x = `Pairwise comparison` , y = `Fold change`, color = Day)) +
geom_point(position = position_dodge(width = 2), size = 3) +
geom_errorbar(
aes(ymin = lower.CL, ymax = upper.CL),
width = .1,
position = position_dodge(width = 2),
size = 1
) +
scale_color_manual(
values = c(
"#ABCFE5",
"#93C4DE",
"#78B5D8",
"#60A6D1",
"#4A97C9",
"#3787C0",
"#2575B7",
"#1664AB",
"#09539D",
"#084185"
)
) +
theme_bw() + theme(panel.grid.major = element_line(colour = "grey100")) +
scale_y_continuous(breaks = log10(c(0.5, 1, 3, 10, 30, 100)),
labels = c(0.5, 1, 3, 10, 30, 100)) +
geom_hline(yintercept = 0,
color = "darkcyan",
lty = "dashed") +
labs(x = "Prime-boost dosing interval",
y = expression(Estimated ~ fold ~ change ~ of ~ S2P ~ IgG ~ titer)) +
facet_grid(`Pairwise comparison` ~ dose,
scales = "free",
space = "free") +
coord_flip() +
theme(
axis.text = element_text(size = 18),
axis.title = element_text(size = 18),
strip.text = element_text(size = 18),
legend.text = element_text(size = 16,
margin = margin(t = 10)),
legend.title = element_text(size = 18),
plot.title = element_text(size = 20),
strip.text.y = element_blank(),
legend.position = "right"
)
Session information
sessionInfo()## R version 4.2.1 (2022-06-23)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] cenGAM_0.5.3 gratia_0.7.3 cowplot_1.1.1 DT_0.24
## [5] pzfx_0.3.0 ggpubr_0.4.0 formattable_0.2.1 lme4_1.1-30
## [9] Matrix_1.5-1 ggeffects_1.1.3 emmeans_1.8.0 mgcv_1.8-40
## [13] nlme_3.1-159 readxl_1.4.1 forcats_0.5.2 stringr_1.4.1
## [17] dplyr_1.0.10 purrr_0.3.4 readr_2.1.2 tidyr_1.2.0
## [21] tibble_3.1.8 ggplot2_3.3.6 tidyverse_1.3.2
##
## loaded via a namespace (and not attached):
## [1] pbkrtest_0.5.1 fs_1.5.2 lubridate_1.8.0
## [4] insight_0.18.2 RColorBrewer_1.1-3 httr_1.4.4
## [7] tools_4.2.1 backports_1.4.1 bslib_0.4.0
## [10] sjlabelled_1.2.0 utf8_1.2.2 R6_2.5.1
## [13] DBI_1.1.3 colorspace_2.0-3 withr_2.5.0
## [16] tidyselect_1.1.2 compiler_4.2.1 cli_3.3.0
## [19] rvest_1.0.3 xml2_1.3.3 labeling_0.4.2
## [22] sass_0.4.2 scales_1.2.1 mvtnorm_1.1-3
## [25] mvnfast_0.2.7 digest_0.6.29 minqa_1.2.4
## [28] rmarkdown_2.16 pkgconfig_2.0.3 htmltools_0.5.3
## [31] highr_0.9 dbplyr_2.2.1 fastmap_1.1.0
## [34] htmlwidgets_1.5.4 rlang_1.0.5 rstudioapi_0.14
## [37] farver_2.1.1 jquerylib_0.1.4 generics_0.1.3
## [40] jsonlite_1.8.0 crosstalk_1.2.0 car_3.1-0
## [43] googlesheets4_1.0.1 magrittr_2.0.3 patchwork_1.1.2
## [46] Rcpp_1.0.9 munsell_0.5.0 fansi_1.0.3
## [49] abind_1.4-5 lifecycle_1.0.1 stringi_1.7.8
## [52] yaml_2.3.5 carData_3.0-5 MASS_7.3-58.1
## [55] grid_4.2.1 parallel_4.2.1 crayon_1.5.1
## [58] lattice_0.20-45 haven_2.5.1 splines_4.2.1
## [61] hms_1.1.2 knitr_1.40 pillar_1.8.1
## [64] boot_1.3-28 estimability_1.4.1 ggsignif_0.6.3
## [67] reprex_2.0.2 glue_1.6.2 evaluate_0.16
## [70] modelr_0.1.9 vctrs_0.4.1 nloptr_2.0.3
## [73] tzdb_0.3.0 cellranger_1.1.0 gtable_0.3.1
## [76] assertthat_0.2.1 cachem_1.0.6 xfun_0.32
## [79] xtable_1.8-4 broom_1.0.1 coda_0.19-4
## [82] rstatix_0.7.0 googledrive_2.0.0 gargle_1.2.0
## [85] writexl_1.4.0 ellipsis_0.3.2